home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 3: Developer Tools / Linux Cubed Series 3 - Developer Tools.iso / devel / lang / lisp / stk-3.002 / stk-3 / STk-3.1 / Mp / fgmp-1.0b5 / gmp.c next >
Encoding:
C/C++ Source or Header  |  1993-12-09  |  37.9 KB  |  1,752 lines
  1. /*
  2.  * FREE GMP - a public domain implementation of a subset of the 
  3.  *           gmp library
  4.  *
  5.  * I hearby place the file in the public domain.
  6.  *
  7.  * Do whatever you want with this code. Change it. Sell it. Claim you
  8.  *  wrote it. 
  9.  * Bugs, complaints, flames, rants: please send email to 
  10.  *    Mark Henderson <markh@wimsey.bc.ca>
  11.  * I'm already aware that fgmp is considerably slower than gmp
  12.  *
  13.  * CREDITS:
  14.  *  Paul Rouse <par@r-cube.demon.co.uk> - generic bug fixes, mpz_sqrt and
  15.  *    mpz_sqrtrem, and modifications to get fgmp to compile on a system 
  16.  *    with int and long of different sizes (specifically MSDOS,286 compiler)
  17.  *  Also see the file "notes" included with the fgmp distribution, for
  18.  *    more credits.
  19.  *
  20.  * VERSION 1.0 - beta 5
  21.  */
  22.  
  23. #include "gmp.h"
  24. #ifndef NULL
  25. #define NULL ((void *)0)
  26. #endif
  27.  
  28. #define iabs(x) ((x>0) ? (x) : (-x))
  29. #define imax(x,y) ((x>y)?x:y)
  30. #define LONGBITS (sizeof(mp_limb) * 8)
  31. #define DIGITBITS (LONGBITS - 2)
  32. #define HALFDIGITBITS ((LONGBITS -2)/2)
  33. #ifndef init
  34. #define init(g) { g = (MP_INT *)malloc(sizeof(MP_INT));  mpz_init(g); }
  35. #endif
  36. #ifndef clear
  37. #define clear(g) { mpz_clear(g); free(g); }
  38. #endif
  39. #ifndef even
  40. #define even(a) (!((a)->p[0] & 1))
  41. #endif
  42. #ifndef odd
  43. #define odd(a) (((a)->p[0] & 1))
  44. #endif
  45.  
  46.  
  47. #ifndef B64
  48. /* 
  49.  * The values below are for 32 bit machines (i.e. machines with a 
  50.  *  32 bit long type)
  51.  * You'll need to change them, if you're using something else
  52.  * If DIGITBITS is odd, see the comment at the top of mpz_sqrtrem
  53.  */
  54. #define LMAX 0x3fffffffL
  55. #define LC   0xc0000000L
  56. #define OVMASK 0x2
  57. #define CMASK (LMAX+1)
  58. #define FC ((double)CMASK)
  59. #define HLMAX 0x7fffL
  60. #define HCMASK (HLMAX + 1)
  61. #define HIGH(x) (((x) & 0x3fff8000L) >> 15)
  62. #define LOW(x)  ((x) & 0x7fffL)
  63. #else
  64.  
  65. /* 64 bit long type */
  66. #define LMAX 0x3fffffffffffffffL
  67. #define LC 0xc000000000000000L
  68. #define OVMASK 0x2
  69. #define CMASK (LMAX+1)
  70. #define FC ((double)CMASK)
  71. #define HLMAX 0x7fffffffL
  72. #define HCMASK (HLMAX + 1)
  73. #define HIGH(x) (((x) & 0x3fffffff80000000L) >> 31)
  74. #define LOW(x) ((x) & 0x7fffffffL)
  75.  
  76. #endif
  77.  
  78. #define hd(x,i)  (((i)>=2*(x->sz))? 0:(((i)%2) ? HIGH((x)->p[(i)/2]) \
  79.     : LOW((x)->p[(i)/2])))
  80. #define dg(x,i) ((i < x->sz) ? (x->p)[i] : 0)
  81.  
  82. #ifdef __GNUC__
  83. #define INLINE inline
  84. #else
  85. #define INLINE
  86. #endif
  87.  
  88. #define lowdigit(x) (((x)->p)[0])
  89.  
  90. struct is {
  91.     mp_limb v;
  92.     struct is *next;
  93. };
  94.  
  95.  
  96. INLINE static void push(i,sp)
  97. mp_limb i;struct is **sp;
  98. {
  99.     struct is *tmp;
  100.     tmp = *sp;
  101.     *sp = (struct is *) malloc(sizeof(struct is));
  102.     (*sp)->v = i;
  103.     (*sp)->next=tmp;
  104. }
  105.  
  106. INLINE static mp_limb pop(sp)
  107. struct is **sp;
  108. {
  109.     struct is *tmp;
  110.     mp_limb i;
  111.     if (!(*sp))
  112.         return (-1);
  113.     tmp = *sp;
  114.     *sp = (*sp)->next;
  115.     i = tmp->v;
  116.     tmp->v = 0;
  117.     free(tmp);
  118.     return i;
  119. }
  120.  
  121. void fatal(s)
  122. char *s;
  123. {
  124.     fprintf(stderr,"%s\n",s);
  125.     exit(123);
  126. }
  127.  
  128. void mpz_init(s)
  129. MP_INT *s;
  130. {
  131.     s->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
  132. #ifdef DEBUG
  133.     printf("malloc: 8 bytes, %08lx\n", (long)(s->p));
  134. #endif
  135.     if (!(s->p))
  136.         fatal("mpz_init: cannot allocate memory");
  137.     (s->p)[0] = 0;
  138.     (s->p)[1] = 0;
  139.     s->sn=0;
  140.     s->sz=2;
  141. }
  142.  
  143. void mpz_init_set(s,t)
  144. MP_INT *s,*t;
  145. {
  146.     int i;
  147.     s->p = (mp_limb *)malloc(sizeof(mp_limb) * t->sz);
  148.     if (!(s->p))
  149.         fatal("mpz_init: cannot allocate memory");
  150.     for (i=0;i < t->sz ; i++)
  151.         (s->p)[i] = (t->p)[i];
  152.  
  153.     s->sn = t->sn;
  154.     s->sz = t->sz;
  155. }
  156.  
  157. void mpz_init_set_ui(s,v)
  158. MP_INT *s;
  159. unsigned long v;
  160. {
  161.     s->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
  162.     if (!(s->p))
  163.         fatal("mpz_init: cannot allocate memory");
  164.     s->p[0] = (v & LMAX);
  165.     s->p[1] = (v & LC) >> DIGITBITS;
  166.     if (v) 
  167.         s->sn = 1;
  168.     else
  169.         s->sn = 0;
  170.     s->sz = 2;
  171. }
  172.  
  173. void mpz_init_set_si(y,v)
  174. MP_INT *y;
  175. long v;
  176. {
  177.     y->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
  178.     if (!(y->p))
  179.         fatal("mpz_init: cannot allocate memory");
  180.     if (v < 0) {
  181.         y->sn = -1;
  182.         y->p[0] = (-v) & LMAX;
  183.         y->p[1] = ((-v) & LC) >> DIGITBITS;
  184.     }
  185.     else if (v > 0) {
  186.         y->sn = 1;
  187.         y->p[0] = v & LMAX;
  188.         y->p[1] = (v & LC) >> DIGITBITS;
  189.     }
  190.     else {
  191.         y->sn=0;
  192.         y->p[0] = 0;
  193.         y->p[1] = 0;
  194.     }
  195.     y -> sz = 2;
  196. }
  197.  
  198.  
  199.  
  200. void mpz_clear(s)
  201. MP_INT *s;
  202. {
  203. #ifdef DEBUG
  204.     printf("free: %08lx\n", (long)(s->p));
  205. #endif
  206.     if (s->p)
  207.         free(s->p);
  208.     s->p=NULL;
  209.     s->sn=0;
  210.     s->sz=0;
  211. }
  212.  
  213. /* number of leading zero bits in digit */
  214. static int lzb(a)
  215. mp_limb a;
  216. {
  217.     mp_limb i; int j;
  218.     j = 0;
  219.     for (i = ((mp_limb)1 << (DIGITBITS-1)); i && !(a&i) ; j++,i>>=1) 
  220.         ;
  221.     return j;
  222. }
  223.  
  224. void _mpz_realloc(x,size)
  225. MP_INT *x; mp_size size;
  226. {
  227.     if (size > 1 && x->sz < size) {
  228.         int i;
  229. #ifdef DEBUG
  230.     printf("realloc %08lx to size = %ld ", (long)(x->p),(long)(size));
  231. #endif
  232.  
  233.         if (x->p) 
  234.             x->p=(mp_limb *)realloc(x->p,size * sizeof(mp_limb));
  235.         else
  236.             x->p=(mp_limb *)malloc(size * sizeof(mp_limb));
  237. #ifdef DEBUG
  238.     printf("becomes %08lx\n", (long)(x->p));
  239. #endif
  240.  
  241.         if (!(x->p))
  242.             fatal("_mpz_realloc: cannot allocate memory");
  243.         for (i=x->sz;i<size;i++)
  244.             (x->p)[i] = 0;
  245.         x->sz = size;
  246.     }
  247. }
  248.  
  249. void dgset(x,i,n)
  250. MP_INT *x; unsigned int i; mp_limb n;
  251. {
  252.     if (n) {
  253.         if (i >= x->sz) {
  254.             _mpz_realloc(x,(mp_size)(i+1));
  255.             x->sz = i+1;
  256.         }
  257.         (x->p)[i] = n;
  258.     }
  259. }
  260.  
  261. static unsigned int digits(x)
  262. MP_INT *x;
  263. {
  264.     int i;
  265.     for (i = (x->sz) - 1; i>=0 && (x->p)[i] == 0 ; i--) 
  266.         ;
  267.     return i+1;
  268. }
  269.         
  270. /* y = x */
  271. void mpz_set(y,x)
  272. MP_INT *y; MP_INT *x; 
  273. {
  274.     unsigned int i,k = x->sz;
  275.     if (y->sz < k) {
  276.         k=digits(x);
  277.         _mpz_realloc(y,(mp_size)k);
  278.     }
  279.     if (y->sz > x->sz) {
  280.         mpz_clear(y); mpz_init(y); _mpz_realloc(y,(mp_size)(x->sz));
  281.     }
  282.  
  283.     for (i=0;i < k; i++)
  284.         (y->p)[i] = (x->p)[i];
  285.  
  286.     for (;i<y->sz;i++)
  287.         (y->p)[i] = 0;
  288.  
  289.     y->sn = x->sn;
  290. }
  291.  
  292. void mpz_set_ui(y,v)
  293. MP_INT *y; unsigned long v; {
  294.     int i;
  295.     for (i=1;i<y->sz;i++)
  296.         y->p[i] = 0;
  297.     y->p[0] = (v & LMAX);
  298.     y->p[1] = (v & LC) >> (DIGITBITS);
  299.     if (v)
  300.         y->sn=1;
  301.     else
  302.         y->sn=0;
  303. }
  304.  
  305. unsigned long mpz_get_ui(y) 
  306. MP_INT *y; {
  307.     return (y->p[0] | (y->p[1] << DIGITBITS));
  308. }
  309.  
  310. long mpz_get_si(y)
  311. MP_INT *y; {
  312.     if (y->sn == 0)
  313.         return 0;
  314.     else
  315.         return (y->sn * (y->p[0] | (y->p[1] & 1) << DIGITBITS));
  316. }
  317.         
  318. void mpz_set_si(y,v)
  319. MP_INT *y; long v; {
  320.     int i;
  321.     for (i=1;i<y->sz;i++)
  322.         y->p[i] = 0;
  323.     if (v < 0) {
  324.         y->sn = -1;
  325.         y->p[0] = (-v) & LMAX;
  326.         y->p[1] = ((-v) & LC) >> DIGITBITS;
  327.     }
  328.     else if (v > 0) {
  329.         y->sn = 1;
  330.         y->p[0] = v & LMAX;
  331.         y->p[1] = (v & LC) >> DIGITBITS;
  332.     }
  333.     else {
  334.         y->sn=0;
  335.         y->p[0] = 0;
  336.         y->p[1] = 0;
  337.     }
  338. }
  339.  
  340. /* z = x + y, without regard for sign */
  341. static void uadd(z,x,y)
  342. MP_INT *z; MP_INT *x; MP_INT *y;
  343. {
  344.     mp_limb c;
  345.     int i;
  346.     MP_INT *t;
  347.  
  348.     if (y->sz < x->sz) {
  349.         t=x; x=y; y=t;
  350.     }
  351.  
  352.     /* now y->sz >= x->sz */
  353.  
  354.     _mpz_realloc(z,(mp_size)((y->sz)+1));
  355.  
  356.     c=0;
  357.     for (i=0;i<x->sz;i++) {
  358.         if (( z->p[i] = y->p[i] + x->p[i] + c ) & CMASK) {
  359.             c=1; 
  360.             (z->p[i]) &=LMAX;
  361.         }
  362.         else 
  363.             c=0;
  364.     }
  365.     for (;i<y->sz;i++) {
  366.         if ((z->p[i] = (y->p[i] + c)) & CMASK)
  367.             z->p[i]=0;
  368.         else
  369.             c=0;
  370.     }
  371.     (z->p)[y->sz]=c;
  372. }
  373.  
  374. /* z = y - x, ignoring sign */
  375. /* precondition: abs(y) >= abs(x) */
  376. static void usub(z,y,x)
  377. MP_INT *z; MP_INT *y; MP_INT *x; 
  378. {
  379.     int i;
  380.     mp_limb b,m;
  381.     _mpz_realloc(z,(mp_size)(y->sz));
  382.     b=0;
  383.     for (i=0;i<y->sz;i++) {
  384.         m=((y->p)[i]-b)-dg(x,i);
  385.         if (m < 0) {
  386.             b = 1;
  387.             m = LMAX + 1 + m;
  388.         }
  389.         else
  390.             b = 0;
  391.         z->p[i] = m;
  392.     }
  393. }
  394.  
  395. /* compare abs(x) and abs(y) */
  396. static int ucmp(y,x)
  397. MP_INT *y;MP_INT *x;
  398. {
  399.     int i;
  400.     for (i=imax(x->sz,y->sz)-1;i>=0;i--) {
  401.         if (dg(y,i) < dg(x,i))
  402.             return (-1);
  403.         else if (dg(y,i) > dg(x,i))
  404.             return 1;
  405.     }
  406.     return 0;
  407. }
  408.  
  409. int mpz_cmp(x,y)
  410. MP_INT *x; MP_INT *y;
  411. {
  412.     int abscmp;
  413.     if (x->sn < 0 && y->sn > 0)
  414.         return (-1);
  415.     if (x->sn > 0 && y->sn < 0)
  416.         return 1;
  417.     abscmp=ucmp(x,y);
  418.     if (x->sn >=0 && y->sn >=0)
  419.         return abscmp;
  420.     if (x->sn <=0 && y->sn <=0)
  421.         return (-abscmp);
  422. }
  423.     
  424. /* is abs(x) == 0 */
  425. static int uzero(x)
  426. MP_INT *x;
  427. {
  428.     int i;
  429.     for (i=0; i < x->sz; i++)
  430.         if ((x->p)[i] != 0)
  431.             return 0;
  432.     return 1;
  433. }
  434.  
  435. static void zero(x)
  436. MP_INT *x;
  437. {
  438.     int i;
  439.     x->sn=0;
  440.     for (i=0;i<x->sz;i++)
  441.         (x->p)[i] = 0;
  442. }
  443.  
  444.  
  445. /* w = u * v */
  446. void mpz_mul(ww,u,v)
  447. MP_INT *ww;MP_INT *u; MP_INT *v; {
  448.     int i,j;
  449.     mp_limb t0,t1,t2,t3;
  450.     mp_limb cc;
  451.     MP_INT *w;
  452.  
  453.     w = (MP_INT *)malloc(sizeof(MP_INT));
  454.     mpz_init(w);
  455.     _mpz_realloc(w,(mp_size)(u->sz + v->sz));
  456.     for (j=0; j < 2*u->sz; j++) {
  457.         cc = (mp_limb)0;
  458.         t3 = hd(u,j);
  459.             for (i=0; i < 2*v->sz; i++) {
  460.                 t0 = t3 * hd(v,i);
  461.                 t1 = HIGH(t0); t0 = LOW(t0);
  462.                 if ((i+j)%2) 
  463.                     t2 = HIGH(w->p[(i+j)/2]);
  464.                 else
  465.                     t2 = LOW(w->p[(i+j)/2]);
  466.                 t2 += cc;
  467.                 if (t2 & HCMASK) {
  468.                     cc = 1; t2&=HLMAX;
  469.                 }
  470.                 else
  471.                     cc = 0;
  472.                 t2 += t0;
  473.                 if (t2 & HCMASK) {
  474.                     cc++ ; t2&=HLMAX;
  475.                 }
  476.                 cc+=t1;
  477.                 if ((i+j)%2)
  478.                     w->p[(i+j)/2] = LOW(w->p[(i+j)/2]) |
  479.                         (t2 << HALFDIGITBITS);
  480.                 else
  481.                     w->p[(i+j)/2] = (HIGH(w->p[(i+j)/2]) << HALFDIGITBITS) |
  482.                         t2;
  483.                     
  484.             }
  485.         if (cc) {
  486.             if ((j+i)%2) 
  487.                 w->p[(i+j)/2] += cc << HALFDIGITBITS;
  488.             else
  489.                 w->p[(i+j)/2] += cc;
  490.         }
  491.     }
  492.     w->sn = (u->sn) * (v->sn);
  493.     if (w != ww) {
  494.         mpz_set(ww,w);
  495.         mpz_clear(w);
  496.         free(w);
  497.     }
  498. }
  499.  
  500. /* n must be < DIGITBITS */
  501. static void urshift(c1,a,n)
  502. MP_INT *c1;MP_INT *a;int n;
  503. {
  504.     mp_limb cc = 0;
  505.     if (n >= DIGITBITS)
  506.         fatal("urshift: n >= DIGITBITS");
  507.     if (n == 0)
  508.         mpz_set(c1,a);
  509.     else {
  510.         MP_INT c; int i;
  511.         mp_limb rm = (((mp_limb)1<<n) - 1);
  512.         mpz_init(&c); _mpz_realloc(&c,(mp_size)(a->sz));
  513.         for (i=a->sz-1;i >= 0; i--) {
  514.             c.p[i] = ((a->p[i] >> n) | cc) & LMAX;
  515.             cc = (a->p[i] & rm) << (DIGITBITS - n);
  516.         }
  517.         mpz_set(c1,&c);
  518.         mpz_clear(&c);
  519.     }
  520. }
  521.     
  522. /* n must be < DIGITBITS */
  523. static void ulshift(c1,a,n)
  524. MP_INT *c1;MP_INT *a;int n;
  525. {
  526.     mp_limb cc = 0;
  527.     if (n >= DIGITBITS)
  528.         fatal("ulshift: n >= DIGITBITS");
  529.     if (n == 0)
  530.         mpz_set(c1,a);
  531.     else {
  532.         MP_INT c; int i;
  533.         mp_limb rm = (((mp_limb)1<<n) - 1) << (DIGITBITS -n);
  534.         mpz_init(&c); _mpz_realloc(&c,(mp_size)(a->sz + 1));
  535.         for (i=0;i < a->sz; i++) {
  536.             c.p[i] = ((a->p[i] << n) | cc) & LMAX;
  537.             cc = (a->p[i] & rm) >> (DIGITBITS -n);
  538.         }
  539.         c.p[i] = cc;
  540.         mpz_set(c1,&c);
  541.         mpz_clear(&c);
  542.     }
  543. }
  544.  
  545. void mpz_div_2exp(z,x,e) 
  546. MP_INT *z; MP_INT *x; unsigned long e;
  547. {
  548.     short sn = x->sn;
  549.     if (e==0)
  550.         mpz_set(z,x);
  551.     else {
  552.         unsigned long i;
  553.         long digs = (e / DIGITBITS);
  554.         unsigned long bs = (e % (DIGITBITS));
  555.         MP_INT y; mpz_init(&y);
  556.         _mpz_realloc(&y,(mp_size)((x->sz) - digs));
  557.         for (i=0; i < (x->sz - digs); i++)
  558.             (y.p)[i] = (x->p)[i+digs];
  559.         if (bs) {
  560.             urshift(z,&y,bs);
  561.         }
  562.         else {
  563.             mpz_set(z,&y);
  564.         }
  565.         if (uzero(z))
  566.             z->sn = 0;
  567.         else
  568.             z->sn = sn;
  569.         mpz_clear(&y);
  570.     }
  571. }
  572.  
  573. void mpz_mul_2exp(z,x,e) 
  574. MP_INT *z; MP_INT *x; unsigned long e;
  575. {
  576.     short sn = x->sn;
  577.     if (e==0)
  578.         mpz_set(z,x);
  579.     else {
  580.         unsigned long i;
  581.         long digs = (e / DIGITBITS);
  582.         unsigned long bs = (e % (DIGITBITS));
  583.         MP_INT y; mpz_init(&y);
  584.         _mpz_realloc(&y,(mp_size)((x->sz)+digs));
  585.         for (i=digs;i<((x->sz) + digs);i++)
  586.             (y.p)[i] = (x->p)[i - digs];
  587.         if (bs) {
  588.             ulshift(z,&y,bs);
  589.         }
  590.         else {
  591.             mpz_set(z,&y);
  592.         }
  593.         z->sn = sn;
  594.         mpz_clear(&y);
  595.     }
  596. }
  597.  
  598. void mpz_mod_2exp(z,x,e) 
  599. MP_INT *z; MP_INT *x; unsigned long e;
  600. {
  601.     short sn = x->sn;
  602.     if (e==0)
  603.         mpz_set_ui(z,0L);
  604.     else {
  605.         unsigned long i;
  606.         MP_INT y;
  607.         long digs = (e / DIGITBITS);
  608.         unsigned long bs = (e % (DIGITBITS));
  609.         if (digs > x->sz || (digs == (x->sz) && bs)) {
  610.             mpz_set(z,x);
  611.             return;
  612.         }
  613.             
  614.         mpz_init(&y);
  615.         if (bs)
  616.             _mpz_realloc(&y,(mp_size)(digs+1));
  617.         else
  618.             _mpz_realloc(&y,(mp_size)digs);
  619.         for (i=0; i<digs; i++)
  620.             (y.p)[i] = (x->p)[i];
  621.         if (bs) {
  622.             (y.p)[digs] = ((x->p)[digs]) & (((mp_limb)1 << bs) - 1);
  623.         }
  624.         mpz_set(z,&y);
  625.         if (uzero(z))
  626.             z->sn = 0;
  627.         else
  628.             z->sn = sn;
  629.         mpz_clear(&y);
  630.     }
  631. }
  632.     
  633.  
  634. /* internal routine to compute x/y and x%y ignoring signs */
  635. static void udiv(qq,rr,xx,yy)
  636. MP_INT *qq; MP_INT *rr; MP_INT *xx; MP_INT *yy;
  637. {
  638.     MP_INT *q, *x, *y, *r;
  639.     int ns,xd,yd,j,f,i,ccc;
  640.     mp_limb zz,z,qhat,b,u,m;
  641.     ccc=0;
  642.  
  643.     if (uzero(yy))
  644.         return;
  645.     q = (MP_INT *)malloc(sizeof(MP_INT));
  646.     r = (MP_INT *)malloc(sizeof(MP_INT));
  647.     x = (MP_INT *)malloc(sizeof(MP_INT)); y = (MP_INT *)malloc(sizeof(MP_INT));
  648.     if (!x || !y || !q || !r)
  649.         fatal("udiv: cannot allocate memory");
  650.     mpz_init(q); mpz_init(x);mpz_init(y);mpz_init(r);
  651.     _mpz_realloc(x,(mp_size)((xx->sz)+1));
  652.     yd = digits(yy);
  653.     ns = lzb(yy->p[yd-1]);
  654.     ulshift(x,xx,ns);
  655.     ulshift(y,yy,ns);
  656.     xd = digits(x); 
  657.     _mpz_realloc(q,(mp_size)xd);
  658.     xd*=2; yd*=2;
  659.     z = hd(y,yd-1);;
  660.     for (j=(xd-yd);j>=0;j--) {
  661. #ifdef DEBUG
  662.     printf("y=");
  663.     for (i=yd-1;i>=0;i--)
  664.         printf(" %04lx", (long)hd(y,i));
  665.     printf("\n");
  666.     printf("x=");
  667.     for (i=xd-1;i>=0;i--)
  668.         printf(" %04lx", (long)hd(x,i));
  669.     printf("\n");
  670.     printf("z=%08lx\n",(long)z);
  671. #endif
  672.         
  673.         if (z == LMAX)
  674.             qhat = hd(x,j+yd);
  675.         else {
  676.             qhat = ((hd(x,j+yd)<< HALFDIGITBITS) + hd(x,j+yd-1)) / (z+1);
  677.         }
  678. #ifdef DEBUG
  679.     printf("qhat=%08lx\n",(long)qhat);
  680.     printf("hd=%04lx/%04lx\n",(long)hd(x,j+yd),(long)hd(x,j+yd-1));
  681. #endif
  682.         b = 0; zz=0;
  683.         if (qhat) {
  684.             for (i=0; i < yd; i++) {
  685.                 zz = qhat * hd(y,i);
  686.                 u = hd(x,i+j);
  687.                 u-=b;
  688.                 if (u<0) {
  689.                     b=1; u+=HLMAX+1;
  690.                 }
  691.                 else
  692.                     b=0;
  693.                 u-=LOW(zz);
  694.                 if (u < 0) {
  695.                     b++;
  696.                     u+=HLMAX+1;
  697.                 }
  698.                 b+=HIGH(zz);
  699.                 if ((i+j)%2) 
  700.                     x->p[(i+j)/2] = LOW(x->p[(i+j)/2]) | (u << HALFDIGITBITS);
  701.                 else
  702.                     x->p[(i+j)/2] = (HIGH(x->p[(i+j)/2]) << HALFDIGITBITS) | u;
  703.             }
  704.             if (b) {
  705.                 if ((j+i)%2) 
  706.                     x->p[(i+j)/2] -= b << HALFDIGITBITS;
  707.                 else
  708.                     x->p[(i+j)/2] -= b;
  709.             }
  710.         }
  711. #ifdef DEBUG
  712.         printf("x after sub=");
  713.         for (i=xd-1;i>=0;i--)
  714.             printf(" %04lx", (long)hd(x,i));
  715.         printf("\n");
  716. #endif
  717.         for(;;zz++) {
  718.             f=1;
  719.             if (!hd(x,j+yd)) {
  720.                 for(i=yd-1; i>=0; i--) {
  721.                     if (hd(x,j+i) > hd(y,i)) {
  722.                         f=1;
  723.                         break;
  724.                     }
  725.                     if (hd(x,j+i) < hd(y,i)) {
  726.                         f=0;
  727.                         break;
  728.                     }
  729.                 }
  730.             }
  731.             if (!f)
  732.                 break;
  733.             qhat++;
  734.             ccc++;
  735. #ifdef DEBUG
  736.             printf("corrected qhat = %08lx\n", (long)qhat);
  737. #endif
  738.             b=0;
  739.             for (i=0;i<yd;i++) {
  740.                 m = hd(x,i+j)-hd(y,i)-b;
  741.                 if (m < 0) {
  742.                     b = 1;
  743.                     m = HLMAX + 1 + m;
  744.                 }
  745.                 else
  746.                     b = 0;
  747.                 if ((i+j)%2) 
  748.                     x->p[(i+j)/2] = LOW(x->p[(i+j)/2]) | (m << HALFDIGITBITS);
  749.                 else
  750.                     x->p[(i+j)/2] = (HIGH(x->p[(i+j)/2]) << HALFDIGITBITS) | m;
  751.             }
  752.             if (b) {
  753.                 if ((j+i)%2) 
  754.                     x->p[(i+j)/2] -= b << HALFDIGITBITS;
  755.                 else
  756.                     x->p[(i+j)/2] -= b;
  757.             }
  758. #ifdef DEBUG
  759.             printf("x after cor=");
  760.             for (i=2*(x->sz)-1;i>=0;i--)
  761.                 printf(" %04lx", (long)hd(x,i));
  762.             printf("\n");
  763. #endif
  764.         }
  765.         if (j%2) 
  766.             q->p[j/2] |= qhat << HALFDIGITBITS;
  767.         else
  768.             q->p[j/2] |= qhat;
  769. #ifdef DEBUG
  770.             printf("x after cor=");
  771.             for (i=xd - 1;i>=0;i--)
  772.                 printf(" %04lx", (long)hd(q,i));
  773.             printf("\n");
  774. #endif
  775.     }
  776.     _mpz_realloc(r,(mp_size)(yy->sz));
  777.     zero(r);
  778.     urshift(r,x,ns);
  779.     mpz_set(rr,r);
  780.     mpz_set(qq,q);
  781.     mpz_clear(x); mpz_clear(y);
  782.     mpz_clear(q);  mpz_clear(r);
  783.     free(q); free(x); free(y); free(r);
  784. }
  785.  
  786. void mpz_add(zz,x,y)
  787. MP_INT *zz;MP_INT *x;MP_INT *y;
  788. {
  789.     int mg;
  790.     MP_INT *z;
  791.     if (x->sn == 0) {
  792.         mpz_set(zz,y);
  793.         return;
  794.     }
  795.     if (y->sn == 0) {
  796.         mpz_set(zz,x);
  797.         return;
  798.     }
  799.     z = (MP_INT *)malloc(sizeof(MP_INT));
  800.     mpz_init(z);
  801.  
  802.     if (x->sn > 0 && y->sn > 0) {
  803.         uadd(z,x,y); z->sn = 1;
  804.     }
  805.     else if (x->sn < 0 && y->sn < 0) {
  806.         uadd(z,x,y); z->sn = -1;
  807.     }
  808.     else {
  809.         /* signs differ */
  810.         if ((mg = ucmp(x,y)) == 0) {
  811.             zero(z);
  812.         }
  813.         else if (mg > 0) {  /* abs(y) < abs(x) */
  814.             usub(z,x,y);    
  815.             z->sn = (x->sn > 0 && y->sn < 0) ? 1 : (-1);
  816.         }
  817.         else { /* abs(y) > abs(x) */
  818.             usub(z,y,x);    
  819.             z->sn = (x->sn < 0 && y->sn > 0) ? 1 : (-1);
  820.         }
  821.     }
  822.     mpz_set(zz,z);
  823.     mpz_clear(z);
  824.     free(z);
  825. }
  826.  
  827. void mpz_add_ui(x,y,n)
  828. MP_INT *x;MP_INT *y; unsigned long n;
  829. {
  830.     MP_INT z;
  831.     mpz_init_set_ui(&z,n);
  832.     mpz_add(x,y,&z);
  833.     mpz_clear(&z);
  834. }
  835.  
  836. int mpz_cmp_ui(x,n)
  837. MP_INT *x; unsigned long n;
  838. {
  839.     MP_INT z; int ret;
  840.     mpz_init_set_ui(&z,n);
  841.     ret=mpz_cmp(x,&z);
  842.     mpz_clear(&z);
  843.     return ret;
  844. }
  845.  
  846. int mpz_cmp_si(x,n)
  847. MP_INT *x; long n;
  848. {
  849.     MP_INT z; int ret;
  850.     mpz_init(&z);
  851.     mpz_set_si(&z,n);
  852.     ret=mpz_cmp(x,&z);
  853.     mpz_clear(&z);
  854.     return ret;
  855. }
  856.  
  857.  
  858. void mpz_mul_ui(x,y,n)
  859. MP_INT *x;MP_INT *y; unsigned long n;
  860. {
  861.     MP_INT z;
  862.     mpz_init_set_ui(&z,n);
  863.     mpz_mul(x,y,&z);
  864.     mpz_clear(&z);
  865. }
  866.     
  867.     
  868. /* z = x - y  -- just use mpz_add - I'm lazy */
  869. void mpz_sub(z,x,y)
  870. MP_INT *z;MP_INT *x; MP_INT *y;
  871. {
  872.     MP_INT u;
  873.     mpz_init(&u);
  874.     mpz_set(&u,y);
  875.     u.sn = -(u.sn);
  876.     mpz_add(z,x,&u);
  877.     mpz_clear(&u);
  878. }
  879.  
  880. void mpz_sub_ui(x,y,n)
  881. MP_INT *x;MP_INT *y; unsigned long n;
  882. {
  883.     MP_INT z;
  884.     mpz_init_set_ui(&z,n);
  885.     mpz_sub(x,y,&z);
  886.     mpz_clear(&z);
  887. }
  888.  
  889. void mpz_div(q,x,y)
  890. MP_INT *q; MP_INT *x; MP_INT *y;
  891. {
  892.     MP_INT r; 
  893.     short sn1 = x->sn, sn2 = y->sn;
  894.     mpz_init(&r);
  895.     udiv(q,&r,x,y);
  896.     q->sn = sn1*sn2;
  897.     if (uzero(q))
  898.         q->sn = 0;
  899.     mpz_clear(&r);
  900. }
  901.  
  902. void mpz_mdiv(q,x,y)
  903. MP_INT *q; MP_INT *x; MP_INT *y;
  904. {
  905.     MP_INT r; 
  906.     short sn1 = x->sn, sn2 = y->sn, qsign;
  907.     mpz_init(&r);
  908.     udiv(q,&r,x,y);
  909.     qsign = q->sn = sn1*sn2;
  910.     if (uzero(q))
  911.         q->sn = 0;
  912.     /* now if r != 0 and q < 0 we need to round q towards -inf */
  913.     if (!uzero(&r) && qsign < 0) 
  914.         mpz_sub_ui(q,q,1L);
  915.     mpz_clear(&r);
  916. }
  917.  
  918. void mpz_mod(r,x,y)
  919. MP_INT *r; MP_INT *x; MP_INT *y;
  920. {
  921.     MP_INT q;
  922.     short sn = x->sn;
  923.     mpz_init(&q);
  924.     if (x->sn == 0) {
  925.         zero(r);
  926.         return;
  927.     }
  928.     udiv(&q,r,x,y);
  929.     r->sn = sn;
  930.     if (uzero(r))
  931.         r->sn = 0;
  932.     mpz_clear(&q);
  933. }
  934.  
  935. void mpz_divmod(q,r,x,y)
  936. MP_INT *q; MP_INT *r; MP_INT *x; MP_INT *y;
  937. {
  938.     short sn1 = x->sn, sn2 = y->sn;
  939.     if (x->sn == 0) {
  940.         zero(q);
  941.         zero(r);
  942.         return;
  943.     }
  944.     udiv(q,r,x,y);
  945.     q->sn = sn1*sn2;
  946.     if (uzero(q))
  947.         q->sn = 0;
  948.     r->sn = sn1;
  949.     if (uzero(r))
  950.         r->sn = 0;
  951. }
  952.  
  953. void mpz_mmod(r,x,y)
  954. MP_INT *r; MP_INT *x; MP_INT *y;
  955. {
  956.     MP_INT q;
  957.     short sn1 = x->sn, sn2 = y->sn;
  958.     mpz_init(&q);
  959.     if (sn1 == 0) {
  960.         zero(r);
  961.         return;
  962.     }
  963.     udiv(&q,r,x,y);
  964.     if (uzero(r)) {
  965.         r->sn = 0;
  966.         return;
  967.     }
  968.     q.sn = sn1*sn2;
  969.     if (q.sn > 0) 
  970.         r->sn = sn1;
  971.     else if (sn1 < 0 && sn2 > 0) {
  972.         r->sn = 1;
  973.         mpz_sub(r,y,r);
  974.     }
  975.     else {
  976.         r->sn = 1;
  977.         mpz_add(r,y,r);
  978.     }
  979. }
  980.  
  981. void mpz_mdivmod(q,r,x,y)
  982. MP_INT *q;MP_INT *r; MP_INT *x; MP_INT *y;
  983. {
  984.     short sn1 = x->sn, sn2 = y->sn, qsign;
  985.     if (sn1 == 0) {
  986.         zero(q);
  987.         zero(r);
  988.         return;
  989.     }
  990.     udiv(q,r,x,y);
  991.     qsign = q->sn = sn1*sn2;
  992.     if (uzero(r)) {
  993.         /* q != 0, since q=r=0 would mean x=0, which was tested above */
  994.         r->sn = 0;
  995.         return;
  996.     }
  997.     if (q->sn > 0) 
  998.         r->sn = sn1;
  999.     else if (sn1 < 0 && sn2 > 0) {
  1000.         r->sn = 1;
  1001.         mpz_sub(r,y,r);
  1002.     }
  1003.     else {
  1004.         r->sn = 1;
  1005.         mpz_add(r,y,r);
  1006.     }
  1007.     if (uzero(q))
  1008.         q->sn = 0;
  1009.     /* now if r != 0 and q < 0 we need to round q towards -inf */
  1010.     if (!uzero(r) && qsign < 0) 
  1011.         mpz_sub_ui(q,q,1L);
  1012. }
  1013.  
  1014. void mpz_mod_ui(x,y,n)
  1015. MP_INT *x;MP_INT *y; unsigned long n;
  1016. {
  1017.     MP_INT z;
  1018.     mpz_init(&z);
  1019.     mpz_set_ui(&z,n);
  1020.     mpz_mod(x,y,&z);
  1021.     mpz_clear(&z);
  1022. }
  1023.  
  1024. void mpz_mmod_ui(x,y,n)
  1025. MP_INT *x;MP_INT *y; unsigned long n;
  1026. {
  1027.     MP_INT z;
  1028.     mpz_init(&z);
  1029.     mpz_set_ui(&z,n);
  1030.     mpz_mmod(x,y,&z);
  1031.     mpz_clear(&z);
  1032. }
  1033.  
  1034. void mpz_div_ui(x,y,n)
  1035. MP_INT *x;MP_INT *y; unsigned long n;
  1036. {
  1037.     MP_INT z;
  1038.     mpz_init_set_ui(&z,n);
  1039.     mpz_div(x,y,&z);
  1040.     mpz_clear(&z);
  1041. }
  1042.  
  1043. void mpz_mdiv_ui(x,y,n)
  1044. MP_INT *x;MP_INT *y; unsigned long n;
  1045. {
  1046.     MP_INT z;
  1047.     mpz_init_set_ui(&z,n);
  1048.     mpz_mdiv(x,y,&z);
  1049.     mpz_clear(&z);
  1050. }
  1051. void mpz_divmod_ui(q,x,y,n)
  1052. MP_INT *q;MP_INT *x;MP_INT *y; unsigned long n;
  1053. {
  1054.     MP_INT z;
  1055.     mpz_init_set_ui(&z,n);
  1056.     mpz_divmod(q,x,y,&z);
  1057.     mpz_clear(&z);
  1058. }
  1059.  
  1060. void mpz_mdivmod_ui(q,x,y,n)
  1061. MP_INT *q;MP_INT *x;MP_INT *y; unsigned long n;
  1062. {
  1063.     MP_INT z;
  1064.     mpz_init_set_ui(&z,n);
  1065.     mpz_mdivmod(q,x,y,&z);
  1066.     mpz_clear(&z);
  1067. }
  1068.  
  1069.  
  1070. /* 2<=base <=36 - this overestimates the optimal value, which is OK */
  1071. unsigned int mpz_sizeinbase(x,base)
  1072. MP_INT *x; int base;
  1073. {
  1074.     int i,j;
  1075.     int bits = digits(x) * DIGITBITS;
  1076.     if (base < 2 || base > 36)
  1077.         fatal("mpz_sizeinbase: invalid base");
  1078.     for (j=0,i=1; i<=base;i*=2,j++)
  1079.         ;
  1080.     return ((bits)/(j-1)+1);
  1081. }
  1082.  
  1083. char *mpz_get_str(s,base,x)
  1084. char *s;  int base; MP_INT *x; {
  1085.     MP_INT xx,q,r,bb;
  1086.     char *p,*t,*ps;
  1087.     int sz = mpz_sizeinbase(x,base);
  1088.     mp_limb d;
  1089.     if (base < 2 || base > 36)
  1090.         return s;
  1091.     t = (char *)malloc(sz+2);
  1092.     if (!t) 
  1093.         fatal("cannot allocate memory in mpz_get_str");
  1094.     if (!s) {
  1095.         s=(char *)malloc(sz+2);
  1096.         if (!s) 
  1097.             fatal("cannot allocate memory in mpz_get_str");
  1098.     }
  1099.     if (uzero(x)) {
  1100.         *s='0';
  1101.         *(s+1)='\0';
  1102.         return s;
  1103.     }
  1104.     mpz_init(&xx); mpz_init(&q); mpz_init(&r); mpz_init(&bb);
  1105.     mpz_set(&xx,x);
  1106.     mpz_set_ui(&bb,(unsigned long)base);
  1107.     ps = s;
  1108.     if (x->sn < 0) {
  1109.         *ps++= '-';
  1110.         xx.sn = 1;
  1111.     }
  1112.     p = t;
  1113.     while (!uzero(&xx)) {
  1114.         udiv(&xx,&r,&xx,&bb);
  1115.         d = r.p[0];
  1116.         if (d < 10) 
  1117.             *p++  = (char) (r.p[0] + '0');
  1118.         else 
  1119.             *p++  = (char) (r.p[0] + -10 + 'a');
  1120.     }
  1121.  
  1122.     p--;
  1123.     for (;p>=t;p--,ps++)
  1124.         *ps = *p;
  1125.     *ps='\0';
  1126.     
  1127.     mpz_clear(&q); mpz_clear(&r); free(t);  
  1128.     return s;
  1129. }
  1130.  
  1131.  
  1132. int mpz_set_str(x,s,base)
  1133. MP_INT *x; char *s; int base;
  1134. {
  1135.     int l,i;
  1136.     int retval = 0;
  1137.     MP_INT t,m,bb;
  1138.     short sn;
  1139.     unsigned int k;
  1140.     mpz_init(&m);
  1141.     mpz_init(&t);
  1142.     mpz_init(&bb);
  1143.     mpz_set_ui(&m,1L);
  1144.     zero(x);
  1145.     while (*s==' ' || *s=='\t' || *s=='\n')
  1146.         s++;
  1147.     if (*s == '-') {
  1148.         sn = -1; s++;
  1149.     }
  1150.     else
  1151.         sn = 1;
  1152.     if (base == 0) {
  1153.         if (*s == '0') {
  1154.             if (*(s+1) == 'x' || *(s+1) == 'X') {
  1155.                 base = 16;
  1156.                 s+=2;   /* toss 0x */
  1157.             }
  1158.             else {
  1159.                 base = 8;
  1160.                 s++;    /* toss 0 */
  1161.             }
  1162.         }
  1163.         else
  1164.             base=10;
  1165.     }
  1166.     if (base < 2 || base > 36)
  1167.         fatal("mpz_set_str: invalid base");
  1168.     mpz_set_ui(&bb,(unsigned long)base);
  1169.     l=strlen(s);
  1170.     for (i = l-1; i>=0; i--) {
  1171.         if (s[i]==' ' || s[i]=='\t' || s[i]=='\n') 
  1172.             continue;
  1173.         if (s[i] >= '0' && s[i] <= '9')
  1174.             k = (unsigned int)s[i] - (unsigned int)'0';
  1175.         else if (s[i] >= 'A' && s[i] <= 'Z')
  1176.             k = (unsigned int)s[i] - (unsigned int)'A'+10;
  1177.         else if (s[i] >= 'a' && s[i] <= 'z')
  1178.             k = (unsigned int)s[i] - (unsigned int)'a'+10;
  1179.         else {
  1180.             retval = (-1);
  1181.             break;
  1182.         }
  1183.         if (k >= base) {
  1184.             retval = (-1);
  1185.             break;
  1186.         }
  1187.         mpz_mul_ui(&t,&m,(unsigned long)k);
  1188.         mpz_add(x,x,&t);
  1189.         mpz_mul(&m,&m,&bb);
  1190. #ifdef DEBUG
  1191.         printf("k=%d\n",k);
  1192.         printf("t=%s\n",mpz_get_str(NULL,10,&t));
  1193.         printf("x=%s\n",mpz_get_str(NULL,10,x));
  1194.         printf("m=%s\n",mpz_get_str(NULL,10,&m));
  1195. #endif
  1196.     }
  1197.     if (x->sn)
  1198.         x->sn = sn;
  1199.     mpz_clear(&m);
  1200.     mpz_clear(&bb);
  1201.     mpz_clear(&t);
  1202.     return retval;
  1203. }
  1204.  
  1205. int mpz_init_set_str(x,s,base)
  1206. MP_INT *x; char *s; int base;
  1207. {
  1208.     mpz_init(x); return mpz_set_str(x,s,base);
  1209. }
  1210.  
  1211. #define mpz_randombyte (rand()& 0xff)
  1212.  
  1213. void mpz_random(x,size)
  1214. MP_INT *x; unsigned int size;
  1215. {
  1216.     unsigned int bits = size * LONGBITS;
  1217.     unsigned int digits = bits/DIGITBITS;
  1218.     unsigned int oflow = bits % DIGITBITS;
  1219.     unsigned int i,j;
  1220.     long t;
  1221.     if (oflow)
  1222.         digits++;
  1223.     _mpz_realloc(x,(mp_size)digits);
  1224.  
  1225.     for (i=0; i<digits; i++) {
  1226.         t = 0;
  1227.         for (j=0; j < DIGITBITS; j+=8)
  1228.             t = (t << 8) | mpz_randombyte; 
  1229.         (x->p)[i] = t & LMAX;
  1230.     }
  1231.     if (oflow) 
  1232.         (x->p)[digits-1] &= (((mp_limb)1 << oflow) - 1);
  1233. }
  1234. void mpz_random2(x,size)
  1235. MP_INT *x; unsigned int size;
  1236. {
  1237.     unsigned int bits = size * LONGBITS;
  1238.     unsigned int digits = bits/DIGITBITS;
  1239.     unsigned int oflow = bits % DIGITBITS;
  1240.     unsigned int i,j;
  1241.     long t;
  1242.     if (oflow)
  1243.         digits++;
  1244.     _mpz_realloc(x,(mp_size)digits);
  1245.  
  1246.     for (i=0; i<digits; i++) {
  1247.         t = 0;
  1248.         for (j=0; j < DIGITBITS; j+=8)
  1249.             t = (t << 8) | mpz_randombyte; 
  1250.         (x->p)[i] = (t & LMAX) % 2;
  1251.     }
  1252.     if (oflow) 
  1253.         (x->p)[digits-1] &= (((mp_limb)1 << oflow) - 1);
  1254. }
  1255.  
  1256. size_t mpz_size(x)
  1257. MP_INT *x;
  1258. {
  1259.     int bits = (x->sz)*DIGITBITS;
  1260.     size_t r;
  1261.     
  1262.     r = bits/LONGBITS;
  1263.     if (bits % LONGBITS)
  1264.         r++;
  1265.     return r;
  1266. }
  1267.  
  1268. void mpz_abs(x,y)
  1269. MP_INT *x; MP_INT *y;
  1270. {
  1271.     if (x!=y)
  1272.         mpz_set(x,y);
  1273.     if (uzero(x))
  1274.         x->sn = 0;
  1275.     else
  1276.         x->sn = 1;
  1277. }
  1278.  
  1279. void mpz_neg(x,y)
  1280. MP_INT *x; MP_INT *y;
  1281. {
  1282.     if (x!=y)
  1283.         mpz_set(x,y);
  1284.     x->sn = -(y->sn);
  1285. }
  1286.  
  1287. void mpz_fac_ui(x,n)
  1288. MP_INT *x; unsigned long n;
  1289. {
  1290.     mpz_set_ui(x,1L);
  1291.     if (n==0 || n == 1)
  1292.         return;
  1293.     for (;n>1;n--)
  1294.         mpz_mul_ui(x,x,n);
  1295. }
  1296.  
  1297. void mpz_gcd(gg,aa,bb)
  1298. MP_INT *gg;
  1299. MP_INT *aa;
  1300. MP_INT *bb;
  1301. {
  1302.     MP_INT *g,*t,*a,*b;
  1303.     int freeg;
  1304.     
  1305.     t = (MP_INT *)malloc(sizeof(MP_INT));
  1306.     g = (MP_INT *)malloc(sizeof(MP_INT));
  1307.     a = (MP_INT *)malloc(sizeof(MP_INT));
  1308.     b = (MP_INT *)malloc(sizeof(MP_INT));
  1309.     if (!a || !b || !g || !t)
  1310.         fatal("mpz_gcd: cannot allocate memory");
  1311.     mpz_init(g); mpz_init(t); mpz_init(a); mpz_init(b);
  1312.     mpz_abs(a,aa); mpz_abs(b,bb); 
  1313.  
  1314.     while (!uzero(b)) {
  1315.         mpz_mod(t,a,b);
  1316.         mpz_set(a,b);
  1317.         mpz_set(b,t);
  1318.     }
  1319.     
  1320.     mpz_set(gg,a);
  1321.     mpz_clear(g); mpz_clear(t); mpz_clear(a); mpz_clear(b);
  1322.     free(g); free(t); free(a); free(b);
  1323. }
  1324.  
  1325. void mpz_gcdext(gg,xx,yy,aa,bb)
  1326. MP_INT *gg,*xx,*yy,*aa,*bb;
  1327. {
  1328.     MP_INT *x, *y, *g, *u, *q;
  1329.  
  1330.     if (uzero(bb)) {
  1331.         mpz_set(gg,aa);
  1332.         mpz_set_ui(xx,1L);
  1333.         if (yy)
  1334.             mpz_set_ui(yy,0L);
  1335.         return;
  1336.     }
  1337.     
  1338.     g = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(g);
  1339.     q = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(q);
  1340.     y = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(y);
  1341.     x = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(x);
  1342.     u = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(u);
  1343.  
  1344.     if (!g || !q || !y || !x || !u)
  1345.         fatal("mpz_gcdext: cannot allocate memory");
  1346.  
  1347.     mpz_divmod(q,u,aa,bb);
  1348.     mpz_gcdext(g,x,y,bb,u);
  1349.     if (yy) {
  1350.         mpz_mul(u,q,y);
  1351.         mpz_sub(yy,x,u);
  1352.     }
  1353.     mpz_set(xx,y);
  1354.     mpz_set(gg,g);
  1355.     mpz_clear(g); mpz_clear(q); mpz_clear(y); mpz_clear(x); mpz_clear(u);
  1356.     free(g); free(q); free(y); free(x); free(u);
  1357. }
  1358.  
  1359.  
  1360. /*
  1361.  *    a is a quadratic residue mod b if
  1362.  *    x^2 = a mod b      has an integer solution x.
  1363.  *
  1364.  *    J(a,b) = if a==1 then 1 else
  1365.  *             if a is even then J(a/2,b) * ((-1)^(b^2-1)/8))
  1366.  *             else J(b mod a, a) * ((-1)^((a-1)*(b-1)/4)))
  1367.  *
  1368.  *  b>0  b odd
  1369.  *
  1370.  *
  1371.  */
  1372.  
  1373. int mpz_jacobi(ac,bc)
  1374. MP_INT *ac, *bc;
  1375. {
  1376.     int sgn = 1;
  1377.     unsigned long c;
  1378.     MP_INT *t,*a,*b; 
  1379.     if (bc->sn <=0)
  1380.         fatal("mpz_jacobi call with b <= 0");
  1381.     if (even(bc))
  1382.         fatal("mpz_jacobi call with b even");
  1383.  
  1384.     init(t); init(a); init(b);
  1385.  
  1386.     /* if ac < 0, then we use the fact that 
  1387.      *  (-1/b)= -1 if b mod 4 == 3
  1388.      *          +1 if b mod 4 == 1
  1389.      */
  1390.  
  1391.     if (mpz_cmp_ui(ac,0L) < 0 && (lowdigit(bc) % 4) == 3)
  1392.         sgn=-sgn;
  1393.  
  1394.     mpz_abs(a,ac); mpz_set(b,bc);
  1395.  
  1396.     /* while (a > 1) { */
  1397.     while(mpz_cmp_ui(a,1L) > 0) {
  1398.  
  1399.         if (even(a)) {
  1400.  
  1401.             /* c = b % 8 */
  1402.             c = lowdigit(b) & 0x07;
  1403.  
  1404.             /* b odd, then (b^2-1)/8 is even iff (b%8 == 3,5) */
  1405.  
  1406.             /* if b % 8 == 3 or 5 */
  1407.             if (c == 3 || c == 5)
  1408.                 sgn = -sgn;
  1409.  
  1410.             /* a = a / 2 */
  1411.             mpz_div_2exp(a,a,1L); 
  1412.  
  1413.         } 
  1414.         else {
  1415.             /* note: (a-1)(b-1)/4 odd iff a mod 4==b mod 4==3 */
  1416.  
  1417.             /* if a mod 4 == 3 and b mod 4 == 3 */
  1418.             if (((lowdigit(a) & 3) == 3) && ((lowdigit(b) & 3) == 3))
  1419.                 sgn = -sgn;
  1420.  
  1421.             /* set (a,b) = (b mod a,a) */
  1422.             mpz_set(t,a); mpz_mmod(a,b,t); mpz_set(b,t);
  1423.         } 
  1424.     }
  1425.  
  1426.     /* if a == 0 then sgn = 0 */
  1427.     if (uzero(a))
  1428.         sgn=0;
  1429.  
  1430.     clear(t); clear(a); clear(b);
  1431.     return (sgn);
  1432. }
  1433.  
  1434. void mpz_and(z,x,y) /* not the most efficient way to do this */
  1435. MP_INT *z,*x,*y;
  1436. {
  1437.     int i,sz;
  1438.     sz = imax(x->sz, y->sz);
  1439.     _mpz_realloc(z,(mp_size)sz);
  1440.     for (i=0; i < sz; i++)
  1441.         (z->p)[i] = dg(x,i) & dg(y,i);
  1442.     if (x->sn < 0 && y->sn < 0)
  1443.         z->sn = (-1);
  1444.     else
  1445.         z->sn = 1;
  1446.     if (uzero(z))
  1447.         z->sn = 0;
  1448. }
  1449.  
  1450. void mpz_or(z,x,y)  /* not the most efficient way to do this */
  1451. MP_INT *z,*x,*y;
  1452. {
  1453.     int i,sz;
  1454.     sz = imax(x->sz, y->sz);
  1455.     _mpz_realloc(z,(mp_size)sz);
  1456.     for (i=0; i < sz; i++)
  1457.         (z->p)[i] = dg(x,i) | dg(y,i);
  1458.     if (x->sn < 0 || y->sn < 0)
  1459.         z->sn = (-1);
  1460.     else
  1461.         z->sn = 1;
  1462.     if (uzero(z))
  1463.         z->sn = 0;
  1464. }
  1465.  
  1466. void mpz_xor(z,x,y)  /* not the most efficient way to do this */
  1467. MP_INT *z,*x,*y;
  1468. {
  1469.     int i,sz;
  1470.     sz = imax(x->sz, y->sz);
  1471.     _mpz_realloc(z,(mp_size)sz);
  1472.     for (i=0; i < sz; i++)
  1473.         (z->p)[i] = dg(x,i) ^ dg(y,i);
  1474.     if ((x->sn <= 0 && y->sn > 0) || (x->sn > 0 && y->sn <=0))
  1475.         z->sn = (-1);
  1476.     else
  1477.         z->sn = 1;
  1478.     if (uzero(z))
  1479.         z->sn = 0;
  1480. }
  1481. void mpz_pow_ui(zz,x,e)
  1482. MP_INT *zz, *x;
  1483. unsigned long e;
  1484. {
  1485.     MP_INT *t;
  1486.     unsigned long mask = (1L<< (LONGBITS-1));
  1487.     
  1488.     if (e==0) {
  1489.         mpz_set_ui(zz,1L);
  1490.         return;
  1491.     }
  1492.  
  1493.     init(t);
  1494.     mpz_set(t,x);
  1495.     for (;!(mask &e); mask>>=1) 
  1496.         ;
  1497.     mask>>=1;
  1498.     for (;mask!=0; mask>>=1) {
  1499.         mpz_mul(t,t,t);
  1500.         if (e & mask)
  1501.             mpz_mul(t,t,x);
  1502.     }
  1503.     mpz_set(zz,t);
  1504.     clear(t);
  1505. }
  1506.  
  1507. void mpz_powm(zz,x,ee,n)
  1508. MP_INT *zz,*x,*ee,*n;
  1509. {
  1510.     MP_INT *t,*e;
  1511.     struct is *stack = NULL;
  1512.     int k,i;
  1513.     
  1514.     if (uzero(ee)) {
  1515.         mpz_set_ui(zz,1L);
  1516.         return;
  1517.     }
  1518.  
  1519.     if (ee->sn < 0) {
  1520.         return;
  1521.     }
  1522.     init(e); init(t); mpz_set(e,ee);
  1523.  
  1524.     for (k=0;!uzero(e);k++,mpz_div_2exp(e,e,1L))
  1525.         push(lowdigit(e) & 1,&stack);
  1526.     k--;
  1527.     i=pop(&stack);
  1528.  
  1529.     mpz_mod(t,x,n);  /* t=x%n */
  1530.  
  1531.     for (i=k-1;i>=0;i--) {
  1532.         mpz_mul(t,t,t); 
  1533.         mpz_mod(t,t,n);  
  1534.         if (pop(&stack)) {
  1535.             mpz_mul(t,t,x); 
  1536.             mpz_mod(t,t,n);
  1537.         }
  1538.     }
  1539.     mpz_set(zz,t);
  1540.     clear(t);
  1541. }
  1542.  
  1543. void mpz_powm_ui(z,x,e,n)
  1544. MP_INT *z,*x,*n; unsigned long e;
  1545. {
  1546.     MP_INT f;
  1547.     mpz_init(&f);
  1548.     mpz_set_ui(&f,e);
  1549.     mpz_powm(z,x,&f,n);
  1550.     mpz_clear(&f);
  1551. }
  1552.     
  1553.     
  1554.  
  1555. /* Miller-Rabin */
  1556. static int witness(x,n)
  1557. MP_INT *x, *n;
  1558. {
  1559.     MP_INT *t,*e, *e1;
  1560.     struct is *stack = NULL;
  1561.     int trivflag = 0;
  1562.     int k,i;
  1563.     
  1564.     init(e1); init(e); init(t); mpz_sub_ui(e,n,1L); mpz_set(e1,e);
  1565.  
  1566.     for (k=0;!uzero(e);k++,mpz_div_2exp(e,e,1L))
  1567.         push(lowdigit(e) & 1,&stack);
  1568.     k--;
  1569.     i=pop(&stack);
  1570.  
  1571.     mpz_mod(t,x,n);  /* t=x%n */
  1572.  
  1573.     for (i=k-1;i>=0;i--) {
  1574.         trivflag = !mpz_cmp_ui(t,1L) || !mpz_cmp(t,e1);
  1575.         mpz_mul(t,t,t); mpz_mod(t,t,n);  
  1576.         if (!trivflag && !mpz_cmp_ui(t,1L)) {
  1577.             clear(t); clear(e); clear(e1);
  1578.             return 1;
  1579.         }
  1580.             
  1581.         if (pop(&stack)) {
  1582.             mpz_mul(t,t,x); 
  1583.             mpz_mod(t,t,n);
  1584.         }
  1585.     }
  1586.     if (mpz_cmp_ui(t,1L))  { /* t!=1 */
  1587.         clear(t); clear(e); clear(e1);
  1588.         return 1;
  1589.     }
  1590.     else {
  1591.         clear(t); clear(e); clear(e1);
  1592.         return 0;
  1593.     }
  1594. }
  1595.  
  1596. unsigned long smallp[] = {2,3,5,7,11,13,17};
  1597. int mpz_probab_prime_p(nn,s)
  1598. MP_INT *nn; int s;
  1599. {   
  1600.     MP_INT *a,*n;
  1601.     int j,k,i;
  1602.     long t;
  1603.  
  1604.     if (uzero(nn))
  1605.         return 0;
  1606.     init(n); mpz_abs(n,nn);
  1607.     if (!mpz_cmp_ui(n,1L)) {
  1608.         clear(n);
  1609.         return 0;
  1610.     }
  1611.     init(a);
  1612.     for (i=0;i<sizeof(smallp)/sizeof(unsigned long); i++) {
  1613.         mpz_mod_ui(a,n,smallp[i]);
  1614.         if (uzero(a)) {
  1615.             clear(a); clear(n); return 0;
  1616.         }
  1617.     }
  1618.     _mpz_realloc(a,(mp_size)(n->sz));
  1619.     for (k=0; k<s; k++) {
  1620.  
  1621.         /* generate a "random" a */
  1622.             for (i=0; i<n->sz; i++) {
  1623.                 t = 0;
  1624.                 for (j=0; j < DIGITBITS; j+=8)
  1625.                     t = (t << 8) | mpz_randombyte; 
  1626.                 (a->p)[i] = t & LMAX;
  1627.             }
  1628.             a->sn = 1;
  1629.             mpz_mod(a,a,n);
  1630.  
  1631.         if (witness(a,n)) {
  1632.             clear(a); clear(n);
  1633.             return 0;
  1634.         }
  1635.     }
  1636.     clear(a);clear(n);
  1637.     return 1;
  1638. }   
  1639.  
  1640.  
  1641. /* Square roots are found by Newton's method, but since we are working
  1642.  * with integers we have to consider carefully the termination conditions.
  1643.  * If we are seeking x = floor(sqrt(a)), the iteration is
  1644.  *     x_{n+1} = floor ((floor (a/x_n) + x_n)/2) == floor ((a + x_n^2)/(2*x))
  1645.  * If eps_n represents the error (exactly, eps_n and sqrt(a) real) in the form:
  1646.  *     x_n = (1 + eps_n) sqrt(a)
  1647.  * then it is easy to show that for a >= 4
  1648.  *     if 0 <= eps_n, then either 0 <= eps_{n+1} <= (eps_n^2)/2
  1649.  *                         or x_{n+1} = floor (sqrt(a))
  1650.  *     if x_n = floor (sqrt(a)), then either x_{n+1} = floor (sqrt(a))
  1651.  *                               or x_{n+1} = floor (sqrt(a)) + 1
  1652.  * Therefore, if the initial approximation x_0 is chosen so that
  1653.  *     0 <= eps_0 <= 1/2
  1654.  * then the error remains postive and monotonically decreasing with
  1655.  *     eps_n <= 2 ^ (-(2^(n+1) - 1))
  1656.  * unless we reach the correct floor(sqrt(a)), after which we may oscillate
  1657.  * between it and the value one higher.
  1658.  * We choose the number of iterations, k, to guarantee
  1659.  *     eps_k sqrt(a) < 1,  so x_k <= floor (sqrt (a)) + 1
  1660.  * so finally x_k is either the required answer or one too big (which we
  1661.  * check by a simple multiplication and comparison).
  1662.  *
  1663.  * The calculation of the initial approximation *assumes* DIGITBITS is even.
  1664.  * It probably always will be, so for now let's leave the code simple,
  1665.  * clear and all reachable in testing!
  1666.  */
  1667. void mpz_sqrtrem (root, rem, a)
  1668.     MP_INT *root;
  1669.     MP_INT *rem;
  1670.     MP_INT *a;
  1671. {
  1672.     MP_INT tmp;
  1673.     MP_INT r;
  1674.     mp_limb z;
  1675.     unsigned long j, h;
  1676.     int k;
  1677.  
  1678.     if (a->sn < 0)
  1679.         /* Negative operand, nothing correct we can do */
  1680.         return;
  1681.  
  1682.     /* Now, a good enough approximation to the root is obtained by writing
  1683.      *     a = z * 2^(2j) + y,  4 <= z <= 15 and 0 <= y < 2^(2j)
  1684.      * then taking
  1685.      *     root = ciel (sqrt(z+1)) * 2^j
  1686.      */
  1687.  
  1688.     for (j = a->sz - 1; j != 0 && (a->p)[j] == 0; j--);
  1689.     z = (a->p)[j];
  1690.     if (z < 4) {
  1691.         if (j == 0) {
  1692.             /* Special case for small operand (main interation not valid) */
  1693.             mpz_set_ui (root, (z>0) ? 1L : 0L);
  1694.             if (rem)
  1695.                 mpz_set_ui (rem, (z>1) ? z-1L : 0L);
  1696.             return;
  1697.         } else {
  1698.             z = (z << 2) + ((a->p)[j-1] >> (DIGITBITS-2));
  1699.             j = j * DIGITBITS - 2;
  1700.         }
  1701.     } else {
  1702.         j = j * DIGITBITS;
  1703.         while (z & ~0xf) {
  1704.             z >>= 2;
  1705.             j += 2;
  1706.         }
  1707.     }
  1708.     j >>= 1;                            /* z and j now as described above */
  1709.     for (k=1, h=4; h < j+3; k++, h*=2); /* 2^(k+1) >= j+3, since a < 2^(2j+4) */
  1710.     mpz_init_set_ui (&r, (z>8) ? 4L : 3L);
  1711.     mpz_mul_2exp (&r, &r, j);
  1712.  
  1713. #ifdef DEBUG
  1714.     printf ("mpz_sqrtrem: z = %lu, j = %lu, k = %d\n", (long)z, j, k);
  1715. #endif
  1716.  
  1717.     /* Main iteration */
  1718.     mpz_init (&tmp);
  1719.     for ( ; k > 0; k--) {
  1720.         mpz_div (&tmp, a, &r);
  1721.         mpz_add (&r, &tmp, &r);
  1722.         mpz_div_2exp (&r, &r, 1L);
  1723.     }
  1724.  
  1725.     /* Adjust result (which may be one too big) and set remainder */
  1726.     mpz_mul (&tmp, &r, &r);
  1727.     if (mpz_cmp (&tmp, a) > 0) {
  1728.         mpz_sub_ui (root, &r, 1L);
  1729.         if (rem) {
  1730.             mpz_sub (rem, a, &tmp);
  1731.             mpz_add (rem, rem, &r);
  1732.             mpz_add (rem, rem, &r);
  1733.             mpz_sub_ui (rem, rem, 1L);
  1734.         }
  1735.     } else {
  1736.         mpz_set (root, &r);
  1737.         if (rem)
  1738.             mpz_sub (rem, a, &tmp);
  1739.     }
  1740.  
  1741.     mpz_clear (&tmp);
  1742.     mpz_clear (&r);
  1743. }
  1744.  
  1745.  
  1746. void mpz_sqrt (root, a)
  1747.     MP_INT *root;
  1748.     MP_INT *a;
  1749. {
  1750.     mpz_sqrtrem (root, NULL, a);
  1751. }
  1752.